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Abstract 

Version 2.0 of CRProp^is public software to model the extra-galactic propa- 
gation of ultra-high energy nuclei of atomic number Z < 26 through structured 
magnetic fields and ambient photon backgrounds taking into account all rele- 
vant particle interactions. CRPropa covers the energy range 6 x 10 16 < E/eV < 
A x 10 22 where A is the nuclear mass number. CRPropa can also be used to 
track secondary 7-rays and neutrinos which allows the study of their link with 
the charged primary nuclei - the so called multi-messenger connection. After a 
general introduction we present several sample applications of current interest 
concerning the physics of extragalactic ultra-high energy radiation. 

Keywords: Ultrahigh energy cosmic rays, Extragalactic magnetic fields. 



1. Introduction and Motivation 

Cosmic rays are ionized atomic nuclei reaching the Earth from outside the 
Solar System with energies that exceed 10 20 eV. Although ultra-high energy 
cosmic rays (UHECRs) were originally discovered in 1939, their sources and 
propagation mechanisms are still a subject of intense research. During the last 
decade significant progress has been made due to the advent of high quality 
and high statistics data from a new generation of large scale observatories. Ob- 
servables of prime interest are the energy spectrum, mass composition and ar- 
rival direction of cosmic rays. A flux-suppression in the energy spectrum above 
E ~ 5 • 10 19 eV has been observed by the HiRes and Pierre Auger Observatories 
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[TJ [2] and possibly also by the Telescope Array |3] indicating either observation 
of the GZK-effect [HE] or the limiting energy of the sources. Moreover, data of 
the Pierre Auger Observatory indicate that the arrival directions of the highest 
energy cosmic rays are correlated with the depth of nearby Active Galactic Nu- 
clei (AGN) or more generally with the nearby extra-galactic matter distribution 

Additionally measurements of the position of the shower maximum and its 
fluctuations by the Pierre Auger Collaboration suggest a significant fraction of 
heavy primaries above 10 19 eV (7]. However, HiRes [8] and preliminary data 
of the Telescope Array [3]), suggest a proton dominance in the same energy 
range. Unfortunately, the limited number of observed events does not yet allow 
the extension of these measurements to the aforementioned cutoff energy. In- 
dependently of the mass composition, it is not uniquely settled yet if this flux 
depression is due to energy loss or maximum energy limitation of the sources. 

Clearly, a better understanding of all of these features and of the effects of 
cosmic ray propagation through the local Universe is mandatory. 

UHECRs do not propagate freely in the Inter Galactic Medium (IGM). Dur- 
ing their propagation they suffer from catastrophic energy losses in reactions 
with the intergalactic background light and are deflected by poorly known mag- 
netic fields. Thus, the effects of propagation alter the cosmic ray spectrum and 
composition injected by sources in the IGM and form the features detected by 
UHECR observatories. In order to establish the origin of UHECRs, it is of 
prime interest to quantitatively understand the imprint of the propagation and 
to disentangle it from the properties of the cosmic rays at their sources. In this 
respect, it is essential to compare the measured UHECR spectrum, composition 
and anisotropy with model predictions. This requires extensive simulations of 
the propagation of UHE nuclei and their secondaries within a given scenario. In 
particular, the observation that UHECRs may consist of a significant fraction of 
heavy nuclei challenges UHECR model predictions and propagation simulations. 
Indeed, compared to the case of ultra-high energy (UHE) nucleons, the propaga- 
tion of nuclei leads to larger deflections in cosmic magnetic fields and additional 
particle interactions have to be taken into account, namely, photodisintegration 
and nuclear decay. 

To provide the community with a versatile simulation tool we present in this 
paper a publicly available Monte Carlo code called CRPropa 2.0 which allows 
one to simulate the propagation of UHE nuclei in realistic one- (ID) and three- 
dimensional (3D) scenarios taking into account all relevant particle interactions 
and magnetic deflections. To this end, we extended the former version 1.4 
of CRPropa, which was restricted to nucleon primaries, to the propagation of 
UHE nuclei. CRPropa 1.4 provided an excellent basis for this effort as many 
of its features could be carried over to the case of UHE nuclei propagation. In 
the present paper, which accompanies the public release of CRPropa 2.0, the 
underlying physical and numerical frameworks of the implementation of nuclei 
propagation are introduced. For technical details the reader is referred to the 
documentation distributed along with this framework. 

This paper is organized as follows: Section [2] starts with a short introduction 
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of the publicly available previous CRPropal.4. The extensions which were 
implemented for nuclei interactions in CRPropa2.0 are the subject of section [3] 
Section [3] describes the general propagation algorithm and in section [5] example 
applications of nuclei propagation with CRPropa are presented. We present a 
short summary and an outlook in Section 6. 

Unless stated otherwise, we use natural units H = c = 1 throughout this 
paper. 

2. Inherited features from CRPropa 1.4 

The previous version 1.4 of CRPropa is a simulation tool aimed at studying 
the propagation of neutrons and protons in the intergalactic medium. It pro- 
vides a one-dimensional (ID) and a three-dimensional (3D) mode. In 3D mode, 
magnetic field- and source distributions can be defined on a 3D grid. This allows 
one to perform simulations in realistic source scenarios with a highly structured 
magnetic field configuration as provided by, e.g., cosmological simulations. In 
ID mode, magnetic fields can be specified as a function of the distance to the 
observer, but their effects are obviously restricted to energy losses of e + e~ pairs 
due to synchrotron radiation within electromagnetic cascades. Furthermore, 
it is possible to specify the cosmological and the source evolution as well as 
the redshift scaling of the background light intensity in ID simulations. All 
important interactions with the cosmic infrared (IRB) and microwave (CMB) 
background light are included, namely, production of electron-positron pairs, 
photopion production and neutron decay. Additionally, CRPropa allows for 
tracking and propagating secondary 7— rays, e + e~ pairs and neutrinos. A mod- 
ule [3] is included that solves the one-dimensional transport equations for elec- 
tromagnetic cascades that are initiated by electrons, positrons or photons taking 
into account single, double and triple pair production as well as up-scattering 
of low energy background photons by inverse Compton scattering. Synchrotron 
radiation along the line of sight can also be simulated. 

Technically, CRPropa is a stand alone object-oriented C++ software pack- 
age. It reads an input file which specifics technical parameters as well as details 
of the simulated "Universe" such as source positions and magnetic fields. The 
CRPropa simulations for a given scenario generate output files of either detected 
events or full UHECR trajectories. 

3. Modeling nuclei interactions in CRPropa 2.0 

Similar to the case of protons, nuclei carry charge and suffer energy losses 
by electron-positron pair production in ambient photon fields. This can occur 
when photon energies boosted into the rest frame of the nucleus are of the 
order of e' ~ 1 MeV. For photon energies at or above the nuclear binding 
energy e' > 8 — 9 MeV, nucleons and light nuclei can be stripped off the nucleus 
(photodisintegration) . Finally, at photon energies exceeding e' ^145 MeV the 
quark structure of free or bound nucleons can be excited to produce mesons 



3 



(photopion production). In these reactions the nucleus can be disrupted and 
unstable elements be produced. Hence, nuclear decay has to be taken into 
account as well. 

In CRPropa2.0 a nucleus with energy E and mass number A is considered a 
superposition of A nucleons with energy E/A. Thus, if one or several nucleons 
are stripped off, the initial energy E will be distributed among the outgoing 
nucleons and nuclei. The ultra-relativistic limit /? — >• 1 (with f3 = v/c of the 
nucleus) is used in CRPropa such that all nuclear products are assigned the 
same velocity vector as the initial particle. This corresponds to a Lorentz factor 
T ~ 10 8 • (i?/10 17 eV) • A^ 1 and a forward collimation within an angle ~ 1/r. 

3.1. Photodisintegration 

Photodisintegration of nuclei has no analogy for free nucleons. Thus, im- 
plementing this new interaction process is mandatory to allow for propagation 
of nuclei within CRPropa. There are many competing photodisintegration pro- 
cesses of different cross sections which need to be accounted for along the path of 
the nucleus in the photon field. Thus, it is important to efficiently describe the 
specific photodisintegration pattern of each propagated nucleus in CRPropa. 

The effects of the propagation of UHE nuclei have first been studied by 
Puget, Stecker and Bredekamp (PSB) [TO]. The approach to model the photo- 
disintegration process chosen in CRPropa 2.0 is similar to what was more re- 
cently discussed in Ref. [TT|. Further details on the photodisintegration within 
CRPropa 2.0 can be found in Ref. [IS]. As target photon fields we shall con- 
sider the CMB and IRB, for which we adopt the more recent parametrization 
developed in [IB] , 

3.1.1. The Photonuclear Cross Sections 

We use the publicly available TALYS framework, version 1.0 [14] to compute 
photodisintegration cross sections. The nuclear models therein are reliable for 
mass numbers A > 12{^] Thus, additional photodisintegration cross sections for 
light nuclei have to added in the modeling. In CRPropa 2.0 TALYS was applied 
to 287 isotopes up to iron (Z = 26)j^] employing nuclear models and settings 
as suggested in Ref. |TT] . The list of isotopes for which the cross sections were 
calculated was generated using data from Ref. [TS]. It is assumed that excited 
nuclei will immediately return to their ground state. Hence, only nuclei in their 
ground states are considered when calculating the cross section^ All cross 
sections were calculated for photon energies IkeV < e' < 250 MeV in the rest 
frame of the nucleus and stored in 500 bins of energy. Knocked out neutrons (n) , 



Nevertheless, TALYS provides results for nuclei with A > 5 and N > 2 where N is the 
number of neutrons. 

3 One could easily extend the framework to nuclei heavier than iron. However, this is not 
done in CRPropa 2.0 because there is no significant contribution of elements with Z > 26 in 
the galactic cosmic ray composition. 

4 In this way we also neglect the photons created in nuclear de-excitations. Note that these 
photons might contribute to the overall flux of photons at energies (E < 0.1 — f EeV). 
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protons (p), deuterium (d), tritium (i), helium-3 ( 3 He) and helium-4 (a) nuclei 
and combinations thereof are considered by TALYS. The corresponding reac- 
tion channels are called exclusive channels. In the mass range of target nuclei 
where TALYS cannot be employed reliably, we use instead other prescriptions, 
as follows: 

• 9 Be, 4 He, 3 He, t and d as given in Ref. |16| . 

• A parametrization of the total photonuclear cross section as function of 
the mass number A is used for 8 Li, 9 Li, 7 Be, 10 Be, n Be, 8 B, 10 B, n B, 
9 C, 10 C, n C as described in |17j . In these cases the loss of one proton 
(neutron) is assumed if the neutron number N < Z (N > Z). For N = Z, 
the loss of one neutron or proton is modeled with equal probability. 

• For 7 Li we use experimental data from Ref. |181 [Hi] and instead of using a 
parametrization we interpolate linearly between the measured data points. 

In total 78449 exclusive channels are taken into account. 

The TALYS output in general agrees reasonably well with available measured 
data and only in rare cases differs up to a factor of 2 for the integrated (total 
absorption) cross sections 

Alternatively, for comparison, the widely-used photodisintegration cross sec- 
tion estimates developed by Puget, Stecker and Bredekamp [TU] can be used. 

In the PSB case a reduced reaction network is implemented involving one 
nucleus for each atomic mass number A up to 56 Fe. Herein, the cross sections 
for one- and two-nucleon dissociation in the photon energy interval e' min < e' < 
30MeV are parametrized by a Gaussian approximation. Different from [10], the 
channel-dependent threshold energies e' min proposed in [5U] are used. In the 
photon energy interval 30 < e' < 150 MeV the cross section is assumed to be 
constant. A comparison of results obtained with the TALYS and PSB cross 
sections is given in section [5. 2| 



3.1.2. Mean Free Path Calculations and Channel Thinning. 

Once the photodisintegration cross section in the nucleus rest frame cr(e') is 



given (cf. Sec. 3.1.1), the energy weighted average cross section 



*(4«) = rAa r^e'a(e')de'. (1) 

\ e max) JO 

is tabulated as a function of e' max = 2Te. 

For a given (7(4,^), for each isotope the mean free path A(r) can be calcu- 
lated as a function of the Lorentz factor according to Ref. [10J , 



A(r)" x = J n(e,z)a( e I nax = 2re)de, (2) 
where n{e, z) is the number density of the isotropic low energy photons per 



energy interval and volume. Variation due to cosmological redshift (cf. Sec. 3.5 1 
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Table 1: Values of the parameters used in Eq. ([TJ and |2j to create the mean free path tables 
for photodisintegration in CRPropa. 



parameter 


CMB 


IRB 


e mi „ (GeV) 


4.00 x 10" 19 




£max (GeV) 


io- n 


10- 7 



is accounted for. For performance reasons, X(E) is tabulated as a function of 
energy. The values of the integration limits e m ; n , e max in Eq. and ([2| are 
listed in Tab. CD 

As the Monte Carlo rapidly slows down with increasing number of exclusive 
channels to be sampled, a thinning procedure was implemented: For each iso- 
tope, we include only the channels with the n largest interaction rates out of the 
N available exclusive channels, such that the sum J^" A~ 1 /A^ > a in at least 
one energy bin. Here, a is the thinning factor and the A" 1 are summed up in 
decreasing order. Furthermore, A^J = A^ 1 is calculated for each isotope. In 
this way a a — 90% channel-thinning reduces the number of photodisintegration 
channels to be tracked from 78449 to 6440. The thinning procedure leads to a 
systematic overestimation of the mean free path in the order of 1% for Lorentz 
factors of the UHECR below 10 10 . This deviation goes up to 10% for Lorentz 
factors above 10 12 (see Fig.[lJ. 

3.2. Photopion Production 

UHE nucleons can produce secondary mesons in interactions with low energy 
photon backgrounds. The most important reaction of this type is the production 
of pions in reactions of UHE protons with the CMB, which leads to the well 
known GZK cut-off [4j [5] at an energy of about Eqzk = 5 • 10 19 eV. 

Nuclei can also produce mesons, albeit with a higher threshold energy of 
Eth ~ A x Eqzk- This is due to the fact that to good approximation the 
center of mass system (CMS) coincides with the rest frame of the nucleus. The 
threshold energy, therefore, does not depend on the total energy of the nucleus 
E, but on the Lorentz factor T cx E/A. Pion production is thus only relevant for 
extremely energetic nuclei. However, it is important to be included to properly 
account for production of secondary UHE photons and neutrinos, as well as for 
the propagation of secondary and primary protons and neutrons. Effectively 
this process leads to an energy scaling of the mean free path Xa,z for photopion 
production of nuclei. The mean free path for pion production on the constituent 
protons A p and neutrons A„ can thus we written as 

X- A ] z {E) = ZX-^{^ +(A-Z)X-i^y (3) 

In CRPropa 2.0 we use Eq. ^ to reduce the mean free path for pion produc- 
tion by nuclei to the one for nucleons which in turn is modeled by the SOPHIA 



G 




Figure 1: The relative deviation of the total mean free path A in photodisintegration reactions 
in the CMB and IRB for the thinned (a = 90%) and un-thinned case is given for all 287 
isotopes (redshift z = 0). 
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package which was already used in CRPropa 1.3. This approximation is suffi- 
cient for our purposes because the pion mass m n rj 140 MeV is much larger than 
the binding energy per nucleon E^/A < 8 MeV so that, above the threshold for 
pion production, the nuclear binding energy can be neglected and the nucleus 
be treated as a collection of free nucleons. Following this argument, we split 
the reaction into four parts, and calculate reactions of protons and neutrons on 
the CMB and on the IRB separately. The nucleus cross section <ja,z for these 
four channels is then given by the cross sections of protons or neutrons times 
the number of the respective nucleons in the nucleus. The mean free path of 
the reaction was then obtained by folding this cross section with the respective 
photon background. 

If a reaction takes place, we treat it as a reaction of a free nucleon. The 
interacting nucleon will suffer energy loss and be stripped off the nucleus. The 
disintegrated nucleon will then be propagated individually, while the produced 
meson will decay leading to secondary leptons, photons or neutrinos which are 
then propagated using the corresponding modules of CRPropa. Both the decay 
of the meson as well as the energy loss of the primary nucleon are calculated by 
using the SOPHIA package [2"T] . 

3.3. Pair production 

Another interaction relevant for the propagation of UHE protons and nuclei 
is the creation of electron positron pairs in the low energy photon backgrounds. 
Both photomeson- and pair-production are less important in terms of energy loss 
of the primary nuclei which is dominated by photodisintegration (c.f. Fig. [3]). 
Pair production is, however, the most important reaction for the creation of 
secondary photons in the TeV range. The mean free path for pair production is 
short, but the energy loss in each individual reaction is small. Thus, we treat 
pair production as a continuous energy loss which for interactions with the CMB 
can be parametrized by |10| 

dE e ^ e ~ 

^f- = 3a cm a T h- 3 Z 2 (m e c 2 k B T)f (T) . (4) 

Here, o~t is the Thomson cross section, m e and m p are the electron and proton 
rest masses, respectively, a om is the fine structure constant, T is the temperature 
of the CMB, and f(T) is a function which depends only on the Lorentz factor 
r and was parametrized by Blumcnthal [22 . One can therefore express the 

energy loss length Ia.z = E {dE A + z /dtj for nuclei in terms of the energy 
loss length for protons l p = E (dEf^ /dt\ , according to 

U,z (T) = ^l p (T) . (5) 

This scaling relation holds for arbitrary target photon backgrounds since the 
prefactor can be traced back to the scaling of the cross section and to the def- 
inition of the energy loss length. Eq. d5l is used in CRPropa 2.0 to generalize 
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the pair production loss rates from protons to nuclei, which in practice is ob- 
tained by integration over the corresponding secondary spectra as parametrized 
by [23]. 

The energy loss for e + e~ pair production is calculated after each timestep At 
and is therefore taken into account at discrete positions and times. CRPropa can 
also propagate secondary electromagnetic cascades initiated by the e + e~ pairs 
or by the 7-rays resulting from 7r°-decay. For the injected secondary spectra 
we use the parametrization given by Kelner and Aharonian [23]. It should be 
noted that, in particular, close to an observer a large time step can degrade 
the accuracy of the propagated spectra, due to the discrete injection of the 
electromagnetic cascade. We refer the reader to the CRPropa 2.0 manual for 
details. 

3-4- Nuclear Decay 

For the propagation of UHECRs, nuclear decay is relevant, if unstable parti- 
cles are produced by photodisintegration or photopion production. On the one 
hand, nuclear decay can change both the nucleus type and its energy, while on 
the other hand it technically ensures that unstable nuclei decay back to stable 
nuclei whose photodisintegration cross sections are known. 

In CRPropa, decays are modeled as a combination of a, f3 ± decays and 
dripping of single nucleons (p, n). The decay length of a nucleus is given by its 
life time r and the Lorentz factor T to be 

Adecay = IV. (6) 

In case of n dripping and a decay, the decay products are assumed to 
inherit the Lorentz factor T from the parent nucleus. This assumption is justified 
since the binding energy per nucleon is small compared to the masses of the 
decay products. The energy of all produced nuclei are, therefore, simply given 
by ' 

Ea\z> = Fm,A',z'- (7) 

In case of (3 ± decay we also use Eq. ([7| and the momenta of e ± and the 
neutrino are calculated from a three body decay (see e.g. |24| ) and are then 
boosted to the simulation frame. 

In CRPRopa the decay channels of the different nuclei as well as their decay 
constants at rest are stored in an internal database. It is based on the NuDat2 
database [TS] and contains 434 different nuclides with mass number A < 56 and 
charge Z < 26. It should be noted that UHECRs, unlike the isotopes in the 
NuDat2 database, are fully ionized. This means that electron capture (EC) is 
not possible for UHECRs and the (3 + decays have to be calculated from the 
EC rates given in the NuDat2 database. Up to the squared matrix elements 
which are the same for EC and /3 + decay, the rates r^c and t7 + for these two 
processes are just proportional to the available phase space of the final state 
products. If Am = itia.z — mA'.z' is the mass difference of the fully ionized 
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Table 2: All isotopes with differing from r by more than 1%. Note that in the case of 36 C7 
the very low branching ratio of 1.6% does not lead to a strong modification of the propagation 
properties. In addition, there are only 3 isotopes where Tp+ is more than 10% larger than t, 
45 Ti, 48 CV and 5i Mn, since 36 CZ and 40 K decay mostly by /3~ decay. Not included in the 
table are 10 isotopes which only decay via electron capture due to their low mass difference 
Am < m e . 



Isotope 




Am/keV 




1.03 


1145 


36 C7 


39.8 


631 




1.56 


993 


43 5c 


1.05 


1710 


45 T j 


1.19 


1551 


47y 


1.01 


2420 


50 V 


1.06 


1696 




1.48 


1144 


49 Cr 


1.03 


2117 


51 Mn 


1.01 


2697 


5i Mn 


4.41 


866 



nuclei and, since the nuclear recoil energy and the electron binding energy can 
be neglected, the kinetic energy of the final state leptons are given by the so- 
called Q-factors Qr+ — Am — m e and Qec = E v = Am + m e = Qp+ + 2m e 
with E v the neutrino energy. Therefore, TjTq oc |V'e(0)| 2 Q| c /7r with |-0 e (O) | 2 = 
(Z/ao) 3 the normalized density of the electron wave function at the nucleus and 
t~1 oc (2/tt 3 ) J^dEf(E) with f(E) = E^E 2 - m^(Am — E) 2 the standard 
j3 + decay differential phase space density per total positron energy E (including 
rest mass) [24]. Thus we get 

tec 2 fa \ 

3 f(E) dE 
Tf}+ ir 2 \z) (Am + m e ) 2 ' 1 ' 

where everything has been expressed in terms of the bare nucleus mass difference 
Am. For all channels involving f3 + decay (including compound channels), we 
compute Tp+ from the lifetime of the dressed nucleus given in NuDat2 r by 
multiplying with Tp+ jr — l+r^+ /tec- The resulting lifetimes Tp+ of all isotopes 
in the database are shown in Fig. [2j In Tab.|2]we list the isotopes for which t^+ 
deviates most from r. 

In Fig. [3] the decay length is shown in comparison to other energy loss pro- 
cesses using the example of 47 Ca. 

3.5. Photon fields and cosmological evolution 

The implementation of photodisintegration and pion production in CR- 
Propa2.0 is based on tabulated mean free path data calculated with the photon 
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Figure 2: Nuclear maps of the internal database. The life time is color-coded, where red 
corresponds to decay times less than 1 s corresponding to a decay length < 0.1 Mpc at a 
Lorentz factor of V = 10 13 . Black denotes effectively stable nuclei with a lifetime larger than 
4 • 10 10 s. This corresponds to a decay length larger than 1/ 'Ho at Lorentz factors of T > 10 7 . 
Nuclei between these two extremes are shown in blue, where darker shades correspond to 
longer decay times. 
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density n(e, z — 0) at a redshift z — (Sec. |3.1.2| . As the photon density n(e, z) 
evolves as a function of z, A = A[r, z] is effectively altered, too. To model this 
change of A as function of z, a scaling function s(z) is used. It approximately re- 
lates ApT, z] at redshift z with the available tabulated data of A[r (1 + z), z = 0] 
at redshift z = 0. For this scaling function s(z), it is assumed that the normal- 
ized spectral shape of the photon field n(e, z) does not change as function of z in 
the comoving cosmological frame. In this approach the evolution of the photon 
number density n(e, z) can be absorbed by a separated evolution factor e(z) 

n(e,z) = {l + zf n^-i-.oj e (z) . (9) 

In the approximation of a redshift independent spectral shape of n(e, z), the 
evolution factor is defined by 

CMB 

, - r r IRB ' ;,0! 

Here, is the intersection energy of the CMB and IRB photon number densities 
n-CMB^i, z) — ^mB(ei) z) in the comoving frame. 

Substitution of Eq. (|9| in Eq. ^ gives the scaling relation for the mean free 
path 

A- J [r,z] = (l + z) 3 e(z) A _1 [r (l + z),z = 0] (11) 

from which one can find s(z) = (1 + z) 3 e(z). It should be noted that this 
result is valid under the assumption that the spectral shape of the IRB does 
not depend on redshift. This approximation does not exactly hold for the IRB, 
due to energy injection in the IGM from galaxy formation. Nevertheless this 
approximation provides a model for the redshift evolution of the IRB which is 
of importance for the production of secondary neutrinos |25| . 




4. Propagation Algorithm and Monte Carlo Approach 

To handle the widely ranging reaction rates of UHE nuclei, a new propaga- 
tion algorithm has been implemented in CRPropa2.0. The main assumption is 
that the mean free paths A are approximately constant during a time step. As 
A = X(E) is in general a function of the UHECR energy E, the numerical step 
size has to be small enough to ensure that no significant energy loss occurs. 

The algorithm works as follows. Given the mean free path Aj for a given in- 
teraction channel of a given nucleus, where i runs over all N possible interaction 
and decay channels, 

1. The inverse total mean free path = YliLi ^ s calculated and a 
distance Axi to the next reaction is selected according to an exponential 
distribution. This is realized by using a uniformly distributed random 
number < r < 1 via 

Axi = -Aln(l-r). (12) 
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Energy (eV) 

Figure 3: Length scales for all interactions of 47 Ca as used by CRPropa. In the plot the dark 
blue (dash-dotted) line denotes the energy loss length of pair production. The red (dotted) 
and green (dashed) lines are the photodisintegration and pionproduction mean free paths 
respectively and the blue (dash-triple-dotted) line shows the decay length. The black (solid) 
line is the total mean free path of the catastrophic energy losses as used in the propagation 
algorithm. 47 Ca was chosen because the half life time of 4.5 days corresponds to a decay 
length comparable to the other interactions. 
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2. The fractional energy loss due to pair production (c.f. Sec. 3.3 1, is limited 
by imposing a maximum step size Ax2 given by 

dx — ^^(E)<5E, (13) 
dx 

where S is the maximal allowed fractional energy loss. 

3. The particle is propagated over a distance 

Ax = min(Axi, Ax 2 , Ax 3 ) , (14) 

where Ax 3 is an upper limit on the propagation step that can be provided 
by the user, with typical values of Ax 3 ^1 — 50 Mpc. This increases 
the accuracy of the calculation of pair production energy losses and the 
accuracy of the secondary pair production spectra. 

4. If Axi = Ax, the particle is propagated over the path length Axi, where 
it performs an interaction. The choice of the specific interaction performed 
by the UHECR is taken by finding the smallest index i for which 

E^r>™, ( 15 ) 

a— 1 

for a uniformly distributed random number < w < 1 . Then the con- 
tinuous energy losses are applied and the algorithm is restarted from the 
new position and the new particles produced in the interaction are added 
to the list of particles to be propagated, leading to a cascade of secondary 
nuclei (see Fig. [4} . A comparison of the exclusive mean free path of the 
different channels A a and the total mean free path for 47 Ca with the pair 
production loss length is shown in Fig. [3] 

If instead Axi > min(Ax2, Ax 3), then continuous energy losses are too 
large to allow for accurate propagation until the next interaction point or 
the user has requested that the maximum step size be smaller than the 
Axi selected in step 1. In this case, the particle is propagated over the 
distance Ax after which continuous losses are applied and the algorithm 
is restarted without performing any interaction. 

From the above description it is clear that if one of the interaction channels 
has a small mean free path X i; the step size of the propagation will adjust itself 
automatically. 

This approach is also applied to select an exclusive channel e.g. in case of 
photodisintegration: If photodisintegration is chosen to be the next reaction by 



the propagation algorithm, the exclusive channel is found by applying Eq. (15 1. 
In analogy, here A to t = (X^ \ rl ) 1 1S the total mean free path for the isotope 
under consideration. The A^ are the mean free path values for the exclusive 
channels of the corresponding isotope. 

If the user chooses to include secondary 7— rays and/or neutrinos, these 
neutral secondaries are propagated over a distance equal to the maximum prop- 
agation distance provided by the user minus the time of their production, such 
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Figure 4: 3D trajectory of an iron nucleus and its hadronic secondaries in the minimum of 
the photodisintegration mean free path (E ~ 1.2 X 10 21 eV) in a high magnetic field region of a 
structured magnetic field (10~ 9 G < B < 10~ 7 G). The iron was injected atx = y = z = 2 Mpc 
along the negative x direction and enters the plotted volume on the top right. Color coded is 
the mass number of the secondary particles. Notice that after photodisintegration the heavy 
nucleus and its secondary particles have the same Lorentz factor T ~ E/A and therefore 
secondary protons are stronger deflected than heavier nuclei due to their higher charge-to- 
mass ratio. 
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that they reach an observer after the maximum propagation time independently 
of the chosen ID or 3D environment. In addition it is possible to inject a mixed 
nuclei composition. Since simple arguments about astrophysical acceleration 
mechanisms indicate that these mixed compositions should be accelerated up 
to a given maximum rigidity R = E/Z in the sources, instead of a maximum 
energy, we included the option to inject up to a given maximum rigidity at the 
source. 

The accuracy of the determination of the arrival direction and arrival time at 
the observer position is of course related to the actual implementation of detec- 
tion and propagation in our algorithm. Besides the numerical error intrinsic to 
the detection algorithm, an additional error is introduced by the choice to take 
continuous losses into account only at the end of the time step and of course by 
the maximum propagation time, which can however be controlled by the user. 
This can be particularly relevant for bursted sources. We refer the interested 
reader to the manual for a deeper discussion on these issues. 

5. Example Applications 

In this section simulations are presented to demonstrate some features of 
CRPropa. All these simulations are restricted to a pure iron or a mixed galactic 
composition injected at the source. For the latter we adopt Ref. [26_, similar to 
the approach of Ref. |27| . The injected power law dN/dE tx E~ a is arbitrarily 
chosen to have a slopaj of a = 2.2. For the 3D simulations, a (75 Mpc) 3 
simulation cube with periodic boundary conditions is defined and filled with 
the large scale structure extragalactic magnetic fields (LSS-EGMF) from the 
cosmological simulations given in Ref. |28| . 

5.1. Completeness of the Photodisintegration Cross Section Tables 

To verify the completeness of the implemented photonuclear cross section 
tables presented in Sec. |3.1.f| i.e. to verify whether cross sections are available for 
all nuclei that occur during propagation, ID simulations have been performed. 
To this purpose, 10 5 iron nuclei were injected with a dN/dE cx E^ 1 spectrum in 
the energy range 1 EeV < E < 56 x 10 3 EeV from a uniform source distribution 
extending up to a distance of 1 Gpc from the observer. In this simulation, all 
particles that were created and propagated within CRPropa were recorded in a 
two dimensional histogram (see Fig. [5]) displaying the frequency of occurrence of 
isotopes in the simulation as a function of their mass and atomic number A, Z. 
In this figure, symbols are given to mark the isotopes for which cross sections 
are available in CRPropa and the type of the marker identifies the source of the 
cross section as listed in Sec. |3.1.l] At the upper right part of Fig. [5]some nuclei 
are created but do not have a photodisintegration cross section assigned. This 



5 The slope a of the source injection spectrum as well as other free parameters of the 
simulations were not optimized in order to reproduce measured UHECR data in these simple 
example applications. 
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Z+0.5 

Figure 5: Isotopes and their frequency of occurrence during the propagation simulated by 
CRPropa for 10 5 injected iron nuclei rectilinearly propagated over a distance up to 1000 Mpc. 
Note that this does not reflect an observable abundance, which would also depend on the time 
of survival of the produced isotope. 
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Figure 6: Average nuclear mass {A) as a function of the distance from the source, resulting 
from the CRPropa photodisintegration cross section tables (see Sec, |3.1,T| l (solid markers) 
and the PSB cross section tables I10I (open markers). The ID simulations assume emission 
of pure iron (black circles) and pure silicon (blue rectangles) with an injection rigidity of 
R = 38.4 EeV. 

(a) Average mass number of the primary cosmic rays, disregarding all nuclei with A < 5. 

(b) Average mass number of all cosmic rays, including secondaries. 



is not problematic since all of them are far off the valley of stability and can 
be expected to decay very quickly. A closer look at the remaining light nuclei 
suggests that five of them cannot be handled due to missing cross section data: 
5 He, 5 Li, 9 B, 7 He and 6 Be. However, all these nuclei have a half-life time smaller 
than 10 -15 s and, hence, are too short-lived to undergo photodisintegration. 
Finally, due to mass and charge loss by pion production, "nuclei" which consist 
only of neutrons or protons can be observed in figure [5] This is a purely technical 
artifact and those nuclei will immediately decay in CRPropa. Thus, this basic 
example simulation suggests that the compilation of photonuclear cross sections 
used in CRPropa 2.0 is complete for an application in UHECR astrophysics. 

5.2. Comparison of the Photodisintegration Cross Section Tables with PSB 
Figure [6] depicts example ID simulations comparing TALYS (Sec. |3.1.T l and 



PSB cross sections [TH]. These simulations show the average mass number (A) 
of UHECRs as a function of the distance from the source for the two cross 
section models under consideration. Here two different cases, namely pure iron 
and pure silicon injection, with the same rigidity at injection R — 38.4 EeV 
(corresponding to injection energies of E =1000 and 538 EeV respectively), are 
considered. Furthermore these two injection cases have been provided both for a 
scenario excluding secondary particles (a) and including secondary particles (b). 
In these simulations, pair production, pion production and redshift evolution 
have been disabled, photodisintegration is considered on the CMB only and the 
cosmic rays are tracked as long as their energy is above 0.1 EeV. In the PSB case, 
decay has been disabled in order to strictly follow the reduced reaction network 
of [in|. As the PSB tables do not provide photodisintegration cross sections 
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for 5 < A < 8, particles that end up in this region are set to immediately 
photodisintegrate to A — 4 plus secondaries. 

In scenario (a), in order to consider only the primary cosmic rays, all nuclei 
with A < 5 have been ignored. This case can be compared with, for instance, 
Fig. 5 of HI]. 

As noted in [IT] , the PSB agreement with experimental data is not as good 
as the one obtained with the Lorentzian parameterization of TALYS. Moreover, 
the PSB case employs a reduced reaction network involving only one nucleus 
for each atomic mass number A up to 56 Fe, whereas 287 nuclides with their 
photodisintegration cross sections have been implemented in CRPropa. Figure[6] 
shows the effect of these differences on the average mass of the primary particle. 
In particular it can be seen that CRPRopa employing TALYS (plus low-mass 
extensions) results on average in a faster photodisintegration rate than PSB 
does. 

In scenario (b) all photodisintegrated particles are included. Due to the 
light secondary particles that are produced with each photodisintegration of the 
primary particle, the average mass number decreases faster in this scenario than 
in scenario (a). This shows the importance of taking secondary particles into 
account when predicting the average mass number. 

A noticeable feature in this scenario is that, for CRPropa, after a certain 
propagation length the average mass number of injected silicon exceeds the 
average mass number of injected iron. This is due to the larger number of light 
secondaries disintegrated off the iron nucleus. This, in combination with the 
cross section dependence on the mass number, can cause a lower total average 
mass at a certain distance, even though the primary cosmic ray still has a higher 
mass on average. 

Furthermore, in this scenario both iron and silicon injection show, at all dis- 
tances, an average mass for the CRPropa tables larger or equal to the average 
mass for the PSB tables. This can be tracked to a difference in the type of 
secondaries that are created. In CRPropa photodisintegration can yield secon- 
daries of mass number up to four (n, p, d, t, 3 He and 4 He). In the PSB case all 
secondaries, with the exception of the reaction 7 + 9 Bc — > 2a + n, have a mass 
number of one, therefore decreasing the average mass number with respect to 
the CRPropa case. 

5.3. ID: Influence of Chemical Composition and Cosmological Evolution 

The advantage of the ID mode in CRPropa is that one can include the 
cosmological as well as the source evolution as function of the distance to the 
observer. We demonstrate this by using two simulations which only differ in 
the composition injected at the sources, namely pure iron injection or a galactic 
mixed composition |26l |2"7] . This allows one to investigate the influence of the 
poorly known initial composition. The parameters of these simulations are as 
follows: UHECRs are injected with an E~ 2 - 2 spectrum up to a rigidity of R = 
384.6 EeV from a continuous source distribution with comoving injection rate 
scaling as (1 + z) 4 up to z max = 2. The cosmological evolution is characterized 
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Figure 7: Results of a ID simulation with CRPropa 2.0 taking into account the cosmological 
expansion and a comoving source evolution scaling as (1 + z) 4 up to z max = 2. 
(Upper panels:) The simulated UHECR flux (black rectangles) has been normalized to the 
Pierre Auger spectrum (red dots) 1291 1301 . The spectra of secondary 7— rays (blue triangles) 
and neutrinos (magenta triangles) have been normalized accordingly. The neutrino flux shown 
is the single-flavor flux, assuming a ratio of 1:1:1. This flux can be compared with the single- 
flavor neutrino limits (black lines) 1311 1321 1331 . Green stars show the isotropic 7-ray flux 
measured by Fermi-LAT [34] . In the left upper panel a galactic mixed composition has been 
injected at the source while in the right upper panel a pure iron composition has been injected. 
(Lower panel:) Abundance of UHE nuclei above 1 EeV after propagation in case of a pure 
iron (blue open circles) and mixed galactic composition (black solid circles) injected. For 
comparison, the original galactic composition (light brown area) is also shown. 
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by a concordance ACDM Universe with a cosmological constant (Q m = 0.3, 
Oa = 0.7) using a Hubble constant of Hq = 72 km s _1 Mpc -1 . 

The results are shown in Fig. [7] The simulated UHECR flux is normalized 
to the spectrum observed by the Pierre Auger Observatory at an energy of 
10 19 EeV [221 ED] and the secondary 7— ray and neutrino fluxes are normalized 
accordingly. The reduction of the flux of observable secondary neutrinos and 
7— rays for pure iron injection shown in Fig. |7(b)] relative to Fig. |7(a) for a mixed 



injection composition is due to photodisintegration dominating with respect to 
photopion production for the heavier composition. The resulting simulated 



abundance of UHE nuclei for E > 1 EeV is shown in Fig 7(c) To compare with 
the propagated composition, the injected galactic composition is also shown. 

Note that photons from nuclear de-excitation during a photodisintegration 
event are currently not taken into account in CRPropa. This might cause a 
moderate increase of the flux of photons at energies below ~ 10 17 eV as discussed 
in e.g. |3S1|3B]. 

5.4- 3D: Continuous Source Distributions following the Large Scale Structure 

A very attractive feature of CRPropa is the possibility to study the effect of 
the presence of a large scale structure extragalactic magnetic field (LSS-EGMF) 
on the UHECR spectrum, composition and anisotropy in 3D simulations. In 
the following example we study how a continuous source distribution following 
the LSS baryon density and deflections in the corresponding LSS-EGMF in the 
scenario of Ref. |28J can influence these quantities. To this end we inject a power 
law of dN/dE oc E^ 1 which we then reweigh to a power law of dN/dE cx E~ 2 - 2 , 
a well known trick to achieve sufficient statistics at high energy. As in the 
previous example, we consider a pure iron and a mixed galactic-like composition 
at injection and for comparison a pure proton composition has been added as 
well. Particles are injected up to a rigidity of R = 384.6 EeV and are tracked as 
long as their energy is above 1 EeV. The detection occurs on a sphere centered 
around the observer, called "sphere around observer" in the code, with a radius 
of ~ 1 Mpc. The observer is placed in a magnetic environment that is similar 
to what is found in the vicinity of our Galaxy. 

Fig. |8(a)] shows that, in case of a pure iron injection, a bump in the UHECR 
spectrum is predicted at ~ 15 EeV, which does not occur if a mixed galactic 
composition or pure proton composition is injected. Furthermore, the simu- 
lated spectra are not strongly affected by the presence of the LSS magnetic field 
or LSS source density. In contrast, the propagated composition of the pure iron 
injection case is affected by deflections, as shown in Fig. |8(b)| since deflections 
increase the propagation path length, thereby enhancing interactions and re- 
ducing the average mass number (A) at detection. Finally, Fig. |8(c) illustrates 



that the pure iron case shows a smaller horizon for energies E ~> 1 EeV when 
compared to the case of an injected galactic or pure proton composition. This is 
mostly due to increased deflections and thus more interactions in case of primary 
iron nuclei. 
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Figure 8: Example simulations of a continuous source distribution which follows the LSS 
density in the scenario of Ref. 28,. Magnetic deflections of the UHECRs in the corresponding 
LSS magnetic fields are taken into account in this 3D simulation. A pure iron composition 
(black rectangles), a galactic composition (red triangles) and a pure proton composition (light 
blue dots) have been injected at the source. 

(a), (b) UHECR flux and average mass number (A) as function of energy. Apart from the three 
cases for injected composition discussed in the text, which are shown including deflections in 
the LSS-EGMF, a pure iron simulation without deflection is shown for comparison (dark blue 
triangles). In the case of the UHECR flux a simulation with pure proton injection without 
deflections and with a flat continuous source distribution (not following the LSS density) has 
been added for comparison as well (magenta stars). The flux is normalized to unity in the 
first bin to allow for a better comparison of the spectral shape. 

(c) Distance of the UHECRs from their source as function of the propagation time for all 
cosmic rays above 1 EeV, including a scenario without deflections (dark blue line). 
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5.5. 3D: Simulating Observables at a given Distance from a Source 

Apart from detecting particles on spheres around the observer, cf. Sec. |5.4| 
CRPropa also allows one to detect particles on spheres around sources. In this 
detection mode a simulated UHECR trajectory crossing a given sphere in a 
3D simulation will be written to the output file. This allows one to study e.g. 
spectrum, composition and anisotropy of the UHECR flux from a given source 
as function of the distance, corresponding to the radii of the spheres. 

As an example a galactic composition is injected from the center of a (75 Mpc) 3 
simulation box filled with the magnetic field configuration as given in the sce- 
nario of Ref. Detection spheres with radii of 4, 8, 16 and 32 Mpc are 
placed around the source. The initial dN/dE cx E~ x spectrum up to a rigidity 
of R = 384.6 EeV is reweighted to a dN/dE oc E~ 2 - 2 spectrum. Only particles 
which have an energy larger than 55 EeV are taken into account. As shown in 
Fig. [9] the flux of cosmic rays at the highest energies is suppressed by particle 
interactions as the distance from the source increases. Furthermore, the cosmic 
ray distribution becomes more anisotropic with increasing distance from the 
source due to increasing deflections in the large scale magnetic field structure. 
This latter effect is exemplified in the sky maps shown in Fig. |9(b)| 



6. Summary and Outlook 

In the present paper we have presented the new version of our UHECR 
propagation code CRPropa, a numerical tool to study the effect of extragalactic 
propagation on the spectrum, chemical composition and distribution of arrival 
directions of UHECRs on Earth. The main new feature introduced in this new 
version 2.0 is the propagation of UHE nuclei, taking also into account their in- 
teractions with the IGM, in particular photodisintegration, which is modeled ac- 
cording to the numerical framework TALYS. As photodisintegration introduced 
many more interaction channels than were present in the previous version of 
CRPropa, we needed to substantially improve the propagation algorithm both 
in efficiency and in accuracy. We also updated the default model for the extra- 
galactic infrared light to a more recent one. CRPropa 2.0 can now be used to 
compute the main observable quantities related to UHECR propagation with 
the accuracy required by present data: particle spectra, mass composition and 
arrival direction on Earth, for highly customizable realizations of the IGM, in- 
cluding source distributions and magnetic fields. In addition, the spectra of 
secondary neutrinos and electromagnetic cascades can be computed down to 
MeV energies. One of the extensions planned for the future is to include a 
module to take into account the effect of galactic magnetic fields on deflections. 
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Figure 9: Spectrum and sky maps of UHECR above 55 EeV after propagating in a (75 Mpc) 3 
simulation cube filled with an LSS-EGMF in the scenario of Ref. 1281 . In this simulation detec- 
tion spheres have been placed around the single source present in the simulated environment 
with the four radii 4, 8, 16, and 32 Mpc and the injection spectrum has been reweighted to a 
dN/dE oc E~ 22 source spectrum. 

(a) Differential spectrum dN/dE multiplied by E 2 as detected at the four spheres. 

(b) Hammer-Aitoff projections of the arrival directions of the simulated UHECR trajectories 
as registered on the different detection spheres around the source. 



24 



by the Deutsche Forschungsgemeinschaft through the collaborative research cen- 
tre SFB 676, by BMBF under grants 05A11GU1 and 05A11PX1, and by the 
"Helmholtz Alliance for Astroparticle Phyics (HAP)" funded by the Initiative 
and Networking Fund of the Helmholtz Association. GS acknowledges sup- 
port from the State of Hamburg, through the Collaborative Research program 
"Connecting Particles with the Cosmos". LM acknowledges support from the 
Alexander von Humboldt foundation. 



References 

[1] The Pierre Auger Collaboration, Observation of the Suppression of the 
Flux of Cosmic Rays above 4xl0 19 eV, Physical Review Letters 101 (2008) 
061101. 

[2] R. U. Abbasi, et al., First Observation of the GZK Cutoff in the HiRcs 
Experiment, Phys. Rev. Lett. 100 (2008) 101101. 

[3] J. N. Matthews, First results from the Telescope Array, 

Nucl.Phys.Proc.Suppl. 212-213 (2011) 79-86. 

[4] K. Greisen, End to the Cosmic-Ray Spectrum?, Phys. Rev. Let. 16 (1966). 

[5] G. T. Zatsepin, V. A. Kuz'min, Upper limit of the spectrum of cosmic rays, 
Sov. Phys. JETP Lett. 4 (1966) 78. Engl. Translation. 

[6] The Pierre Auger Collaboration, Correlation of the Highest-Energy 
Cosmic-Rays with Nearby Extragalactic Objects, Science 318 (2007). 

[7] J. Abraham, et al., Measurement of the Depth of Maximum of Extensive 
Air Showers above 10 18 eV, Phys.Rev.Lett. 104 (2010) 091101. 

[8] R. U. Abbasi, et al., Indications of Proton-Dominated Cosmic-Ray Com- 
position above 1.6 EeV, Physical Review Letters 104 (2010) 161101. 

[9] S. Lee, On the propagation of extragalactic high-energy cosmic and gamma- 
rays, Phys.Rev. D58 (1998) 043004. 

[10] J. Puget, F. Stecker, J. Bredekamp, Photonuclear Interactions of Ultrahigh- 
Energy Cosmic Rays and their Astrophysical Consequences, Astrophys.J. 
205 (1976) 638-654. 

[11] E. Khan, et al., Photodisintegration of ultra-high-energy cosmic rays re- 
visited, Astroparticle Physics 23 (2005) 191-201. 

[12] N. Nierstcnhoefer, On the Origin and Propagation of Ultra-High Energy 
Cosmic Rays, Ph.D. thesis, University of Wuppertal, 2011. 

[13] T. Kneiske, private communication, 2011. University of Hamburg. 



25 



[14] A. Koning, et al., TALYS-1.0, in: Proceedings of the International Con- 
ference on Nuclear Data for Science and Technology, EDP Sciences, 2008, 
pp. 211-214. 

[15] National Nuclear Data Center, NuDat 2 database, See 

http://www.nndc.bnl.gov/nudat2/, 2011. 

[16] J. Rachen, Interaction Processes and Statistical Properties of the Propaga- 
tion of Cosmic Rays in Photon Backgrounds, Ph.D. thesis, Universitat zu 
Bonn, 1996. 

[17] M. V. Kossov, Approximation of photonuclear interaction cross-sections, 
European Physical Journal A 14 (2002) 377-392. 

[18] V. V. Varlamov, et al., PHOTONUCLEAR DATA. PHOTODISIN- 
TEGRATION OF LITHIUM. EVALUATED CROSS SECTIONS OF 
CHANNELS AND REACTIONS., Fotojad.Dannye - Photodisint.of 
Li.Suppl., Moscow 1986, 1986. 

[19] L. A. Kulchitskii, Y. M. Volkov, V. P. Denisov, V. I. Ogurtsov, Energy 
levels of li 7 observed in its photoemission, Izv. Rossiiskoi Akademii Nauk, 
Ser.Fiz. 27 (1963) 1412. 

[20] F. Stecker, M. Salamon, Photodisintegration of ultrahigh-energy cosmic 
rays: A New determination, Astrophys.J. 512 (1999) 521-526. 

[21] A. Miicke, R. Engel, J. P. Rachen, R. J. Protheroe, T. Stanev, Monte Carlo 
simulations of photohadronic processes in astrophysics, Computer Physics 
Communications 124 (2000) 290-314. 

[22] G. R. Blumenthal, Energy loss of high-energy cosmic rays in pair-producing 
collisions with ambient photons, Phys. Rev. D 1 (1970) 1596-1602. 

[23] S. R. Kelner, F. A. Aharonian, Energy spectra of gamma rays, electrons, 
and neutrinos produced at interactions of relativistic protons with low en- 
ergy radiation, prd 78 (2008) 034013. 

[24] J. Basdevant, J. Rich, M. Spiro, Fundamentals in nuclear physics: From 
nuclear structure to cosmology, 2005. 

[25] K. Kotera, D. Allard, A. Olinto, Cosmogenic Neutrinos: parameter space 
and detectabilty from PeV to ZeV, JCAP 1010 (2010) 013. 25 pages, 12 
figures, version to appear in JCAP, minor changes. 

[26] M. A. Duvernois, M. R. Thayer, The Elemental Composition of the Galactic 
Cosmic-Ray Source: ULYSSES High-Energy Telescope Results, Astrophys. 
J. 465 (1996) 982. 

[27] D. Allard, E. Parizot, E. Khan, S. Goriely, A. Olinto, UHE nuclei prop- 
agation and the interpretation of the ankle in the cosmic-ray spectrum, 
Astron.Astrophys. 443 (2005) L29-L32. 



26 



[28] G. Sigl, F. Miniati, T. A. Enfilin, Cosmic Magnetic Fields and Their Influ- 
ence on Ultra-High Energy Cosmic Ray Propagation, Nuclear Physics B 
Proceedings Supplements 136 (2004) 224-233. 

[29] J. Abraham, et al., Measurement of the energy spectrum of cosmic rays 
above 10 18 eV using the Pierre Auger Observatory, Phys.Lett. B685 (2010) 
239-246. 

[30] P. Abreu, et al., The Pierre Auger Observatory I: The Cosmic Ray En- 
ergy Spectrum and Related Measurements (Contribution to the 32 nd ICRC, 
Beijing, arXiv: 1107.4809 [astro-ph.HE]). 

[31] R. Abbasi, et al., Constraints on the Extremely-high Energy Cosmic Neu- 
trino Flux with the IceCube 2008-2009 Data, Phys.Rev. D83 (2011) 092003. 

[32] P. Abreu, et al., A Search for Ultra-High Energy Neutrinos in Highly 
Inclined Events at the Pierre Auger Observatory, Phys.Rev. D84 (2011) 
122005. 

[33] P. Gorham, et al., Observational Constraints on the Ultra-high Energy 
Cosmic Neutrino Flux from the Second Flight of the ANITA Experiment, 
Phys.Rev. D85 (2012) 049901. 

[34] A. Abdo, et al., The Spectrum of the Isotropic Diffuse Gamma-Ray 
Emission Derived From First- Year Fermi Large Area Telescope Data, 
Phys.Rev.Lett. 104 (2010) 101101. 

[35] L. A. Anchordoqui, J. F. Bcacom, H. Goldberg, S. Palomares-Ruiz, 
T. J. Weiler, TeV gamma-rays from photo-disintegration/de-excitation of 
cosmic-ray nuclei, Phys.Rev.Lett. 98 (2007) 121101. 

[36] K. Murase, J. F. Beacom, Very-High-Energy Gamma-Ray Signal from Nu- 
clear Photodisintegration as a Probe of Extragalactic Sources of Ultrahigh- 
Energy Nuclei, Phys.Rev. D82 (2010) 043008. 



27 



